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1 . Introduction 

We consider here finite algorithms for solving the nonlinear program PI : 

minimize f(yi,...,yn) 



m 

subject to X E(iJ)xij=yj; l<j<n, 
i=l 
n 

X Xij=bi; l<i<m, 

j=> 

xjj>0; l<i<m, l<j<n. 



It is assumed in PI that bi>0 for all i and that E(i ,j)>0 for all (i,j), with at least one positive 
E(i,j) in each row i. It is further assumed that f(y) is continuous and decreasing in each of 
its arguments on some convex set S that includes the nonnegative orthant of R". y 
represents the vector (yp. . .,yn)» and similarly x will represent the collection of all mn of 
the xij, etc. PI is a special case of the nonlinear network problem, a problem that is already 
known to possess sufficient special structure to make the development of special purpose 
software attractive(Ahlfeld, Dembo, Mulvey, and Zenios(1987)). Our intention is to 
specialize PI even further by making restrictive assumptions about f(y). The components 
of y will be called "potentials", one for each column of the matrix (E(i,j)). Potentials are 
not logically necessary in defining PI, since the expressions defining y could simply be 
substituted into f(y), but the generous notation will prove useful in the sequel. The 
components of x will be called "allocations." 

PI could be interpreted as a Search Theory problem by assuming that search is to be 
conducted over a fixed time period by m distinct tjqDes of searcher, with bj units of type i 
effort available over the period. Assume that a single target is located with probability pj in 
one of n regions, that region J has area Aj, and that the reward for finding the target in 
region j is Vj Assume funher that the amount of area swept by a searcher of type i in 
region j is Sy. If Xy is the allocation of search effort of type i to region j, then the total area 

in 

swept in region j is X ^^e search in each area is "random" (Koopman(1956)), 

i=l 

then the number of times a target in region j is detected is a Poisson random variable with 






mean 



s 



XijSij 



v= 



/Aj =yj (so let E(i,j)=Sy/Aj), and the probability of no detections is 



n 

exp(-yj). The average "regret" in terms of value not found is then f(y)=X '^jPj^xp(“yj)- 

J=.’ . 

PI is then the problem of allocating m resources to n regions in order to minimize the 
average undetected value. PI can be thought of as a generalization from the case m=l, 
which already has specially tailored solution procedures (Washbum(1981)). 

The US Air Force's Heavy Attack model (Clasen, Graves, and Lu(1974)) is also 
closely related to PI , with resources being aircraft sorties, areas being classes of targets. 
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and potentials being "average number of targets killed were it not for the fact that some of 
the potential will be 'wasted' in attacking targets already dead". In general, PI is applicable 
when n projects must be simultaneously undertaken by employing m resources, y is a 
vector measure of progress, and f(y) is some global, scalar measure of work not 
completed. 

2. Basic Feasible Solutions 

PI has mn potentially positive allocations, but optimal solutions invariably have most 
of these being 0. The reason for this is that optimal solutions need never include cycles of 
positive allocations. 

Definition: If a feasible allocation for PI has an alternating sequence 
jl>il J242»---qLdL+l of columns and rows all distinct except that jL+i=ji, and if ^>0 for 
every consecutive pair in the sequence ((i,j)=(ikjk) or (i,j)=(ik, jk+l) for some k^), then 
the sequence is a cycle. 

Definition: A feasible allocation for PI is conservative if Xij=0 whenever E(i,j)=0. 

Theorem 1: Given the assumptions of section 1, PI has an optimal, conservative 
solution with no cycles. 

Proof: PI has an optimal solution because the objective function regarded as a function 
of X is continuous on a compact set (Bazaraa and Shetty(1979)). An optimal solution x 
can easily be converted to a conservative optimal solution by shifting any offending 
allocation xu from column j to some other column k where E(i,k)>0 (by assumption at 
least one such column k always exists for every i). Assume therefore that x is optimal and 
conservative, but that a cycle exists. In that case construct a modified solution x' by adding 

5(i,j) to Xjj for each consecutive pair (i,j) in the sequence, where 5(i,j) is defined as 
follows: 

5(ii.ji)=z (an arbitrary real number) 

5(ikJk+l)=- 5(ik,jk); 

5(ik+lJk+l)=- 5(ik,jk+l)E(ikJk+l)/E(ik+l jk+l); l^^-l. 

n 

The numbers 5(i,j) are constructed so that ^ x'ij=bi for all i as before. Also 

j=i 

m 

^ x'ijE(iJ)=yj as before, except that the latter sum for j=ji is increased by 
i=l 

5(ii,ji)E(ii,ji)-6(iLji)E(iLji). This difference is proportional to z. If the proportionality 
constant is nonnegative, increase z until x'jj=0 for some (i,j) in the cycle; there will be 

some such (i,j) because the numbers 5ij are proportional to z with nonzero proportionality 
constants that alternate in sign. Otherwise, decrease z to achieve the same end Since f() is 
decreasing in each variable, the x' solution is at least as good as the x solution, with one 
less positive variable. By repeating this operation if necessary, an optimal solution with no 
cycles can be found. QED 

Since an optimal conservative solution with no cycles exists, it makes sense to search 
for the optimal solution amongst solutions with that property. A major advantage of such 
solutions is that there can be at most n+m-1 positive variables xjj. To see this, consider the 
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nonoriented bipartite graph G formed by connecting m "row nodes" to n "column nodes", 
with i being connected to j if and only if xij>0. Such a graph must consist of a number of 
connected components Ti,...Tk, each of which is not connected to any of the others. 

Some of these components may contain a single node and no edges, but every node is in 
exactly one component. Since each component is connected and has no cycles, it is a tree 
(it is convenient to refer to even singleton components as "trees"), so G is a "forest". If 
has rk nodes, it must have exactly rj^l edges (Berge(1962)), even if rk=l. Therefore the 
total number of edges is n+m-K, which is at most n+m-1. 

The graph G described above will be called the "basis" of any conservative solution 
with no cycles. More generally. 

Definition: Any nonoriented bipartite graph G with the following properties is a basis: 

1) The nodes of G are all of the m+n row and column nodes. 

2) All row nodes I for which bi>0 have at least one incident edge. 

3) If E(i,j)=0, then there is no edge (i,j). 

4) There are no cycles. 

The requirement that a solution should have no cycles is sufficient to establish that the 
corresponding basis is a forest, but it is not sufficient to establish a one-to-one relationship 
between bases and solutions. In order to establish such a relationship, further 
assumptions about f() are needed. 

Definition: f() is good if it has the following properties: 

1) f() is differentiable and decreasing on the interior of the convex set S described 
earlier. Let f() be the gradient of f(). 

2) There exists a unique function g(|a) taking values in S such that f (g(ia))+|i=0 for 
all |i>0, where by |i>0 we mean ppO for j=l,...,n. 

3) If C is a nonempty subset of { l,...,n), if |i>0, and if scalar p>0, then the equation 
^|ijgj(z|i)=p has a unique positive solution z. Furthermore z is a continuous function of 



In practice the methods introduced below will be attractive only if the solution z 
required in the last part is easily computed, since the equation will have to be solved many 
times. Here are two examples of good functions: 



jeC 

p for p>0. 



n 



Example 1: f(y)=^ Vjexp(-yj), with S=R". This is the Search Theory function 
j=i 

introduced earlier, except that pj has been incorporated into Vj. In this case gj(|i)=ln(Vj/|ij) 





j=i 



and Vz= VjUj / p+Xl^j ’ defined for p>0. 



jeC J L jeC 



The Kiihn-Tucker conditions for PI are that there should exist and y 



such that 
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KTl) ?.i>^jE(i,j) for all i,j 
KT2) xjj>0 for all i,j 

KT3) ^i=[ajE(i,j) when xij?^0 

m 

KT4) ^ XijE(i,j)=yj for all j 
i=i 
n 

KT5) ^ xij=bj for all i, and 
j=i 

KT6) y=g(^). 

The object now is to use the equality parts of the KT conditions (KT3-KT6) to set up a 
one-to-one correspondence between bases and solutions for good objective functions, after 

which solutions corresponding to bases will be called "basic." The vectors X and )i will 
each be called "multipliers." Basic solutions will have xjj=0 for nonbasic edges (i,j), but 
the following theorem permits one exception to support the pivoting ideas of section 4. 

Theorem 2: For any tree T let R(T) and C(T) be the row nodes and column nodes of T, 

respectively. Let G be a basis composed of trees (Ti,...,Tk), and suppose that ?ii=|ijE(i,j) 
for all edges (i,j) of G. Then 

a) If (x,y) solve KT4-KT5, and if all nonbasic allocations except for xjj are zero, then 
for any component T of G, 



Z Pjyj 


-pjxijE(I,J)5(C(T),J) = 


Z^ibi 


-?iIxij5(R(T),I), 


(1) 


J£C(T) 




_ i£R(T) 







where now 5(W,w) is 1 if w is in the set W, otherwise 0. Furthermore, 

b) If (^,|i,y) is such that (1) holds for every component T of G, then KT4-KT5 have 
a unique solution x for which all nonbasic allocations except xu are zero. 

Proof of part a) Multiply both sides of KT4 by [ij and sum to obtain 

m 

Z Z FjXijE(i,j)= X FjYj- (2) 

j£C(T) i=l j£C(T) 

Since xij=0 when jeC(T) unless either i=I or ieR(T), 

Z Z l^jXijE(i,j) + ^ijxijE(I,J)5(c(T),J) = ^ pjyj. (3) 

jeC(T) i£R(T) j£C(T) 

But |ijE(i,j)=?Li when jeC(T) and ieR(T) by assumption, so 

Z^i Z ^ij + W^IJE(EJ)5 (c(T),J) = ^ pjyj. (4) 

i£R(T) j£C(T) jeCCT) 

Since KT5 holds, 

Z >^ij =bi - xu5({i),D for all ieR(T). (5) 

j£C(T) 

Part a) follows upon substituting (5) into (4). 
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Proof of part b): To avoid unnecessary complication assume xu=0; otherwise bj and yj 
can be adjusted by subtracting xjj and xijE(l,J), respectively. Consider any component T 
of G. We must show that there is exactly one way to assign allocations xjj to the edges of 
T that is consistent with KT4-KT5. This is trivial if T is a singleton; otherwise, since T is 
a tree, there must be at least 1 pendant (connected to exactly one other node) node. If this 
pendant node is row P, let row P be connected to column Q. Then xpQ is the only nonzero 
allocation in row P, and is therefore necessarily bp by KT5. Define y' by 
y’Q=yQ- bpE(P,Q), yj=yj for J?>iQ, and let T' be the tree remaining after node P and edge 

(P,Q) are deleted. Since ?ip=|iQE(P,Q), |iQy'Q=(tQyQ-?ipbp. Therefore 

X ^jy 'i “ S ^ pendant node is column Q, let it be connected to row P. 
jeC(T') ieR(T') 

Since xpQ is the only positive allocation in column Q, necessarily xpQ is yQ/E(P,Q) 
according to KT4. xpQ is well defined because E(P,Q)>0. Let b'p=bp-xpQ, otherwise let 
b’i=bj for i?!:P, and let T' be the tree resulting from deleting column Q and edge (P,Q). 

Since Xp=|IqE(P,Q), Xpb'p=Xpbp-|J.QyQ. Therefore pjyj = ^ ^tjb'i. Since in 

jeC(T') ieR(T') 

either case T' has one less node than T, and since the proposition is true for graphs with 
one node, the theorem is proved by induction. QED 

Multipliers for which Xi=pjE(i,J) on basic edges are easy to detemiine if there is no 
requirement for equation ( 1) to hold. The reason is that if or pj is determined at any 

node, then the same thing is true of all neighbors of the node. Thus X and p can be 
determined in any component T of G by assigning an arbitrary positive value to any node. 

Let (^*,p*) be such multipliers. Then any other set (I^,p) with the same property is some 
scalar multiple of (?i*,p*). Let z>0 be the multiple. Then z must be detennined so that (1) 
holds. After substituting gj(zp*) for yj, (1) is; 

X -p*jxijE(I,J)6(C(T),J)= X - ^i^^'ixijSCRCT),!). (6) 

jeC(T) >eR(T) 

As long as xu<bi, the right hand side of (6) is nonnegative. Equation (6) will then have a 
unique positive solution z as long as f() is good. Given z, KT3-KT6 have the unique 

solution (X,p)=z(?i*,p*), y=g(p), and x the unique solution of KT4—KT5. The x part is 
the desired basic solution corresponding to the basis G-hereafter the G-basic solution. 

Thus the unique G-basic solution can be easily determined as long as xu^bi . 

A G-basic solution need not satisfy KT2, but if it does then it will be called a basic 
feasible solution and G will be called "feasible". Some basic feasible solution is optimal, 
so PI could be solved in principle by examining all bases. In analogy with the Simplex 
Method, however, the object is to avoid this by utilizing a procedure that is directional in 
the sense of considering a sequence of feasible bases, each of which is better than the last. 
Two of these will be considered below. In each case we will prove only that the objective 
function sequence is nonincreasing, rather than decreasing, with equahty being possible in 
degenerate cases. Rather than deal with the distracting issue of degeneracy here, we simply 
assume nondegeneracy for the next two sections. There is every reason to expect that the 
same techniques that handle degeneracy in the Simplex Method will also work here. 
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3 . The Manifold Suboptimization Method 

In this section nonbasic variables will always be 0. 

A G-basic solution is the optimal solution of P2(G): 

mimimize f(y) 

m 

subject to X 
i=l 
n 

xjj<bi; l^<m, 
j=i 

xij=0 unless (i,j)eG. 

To prove this note that the G-basic solution solves KT3-KT6. But KT3-KT6 are the 
Kuhn Tucker conditions for P2(G), and the Kuhn Tucker conditions are sufficient for a 
global minimum because f() is convex as a function of x. P2(G) does not require 
allocations to be nonnegative, so the feasible space is not compact. Absent the constraints 
forcing nonbasic variables to be 0, P2(G) might lack a (finite) optimal solution. 
Nonetheless, the G-basic solution is optimal. 

Given a basis G and a corresponding solution x feasible in PI, suppose KTl is not 
satisfied for some (I,J). The natural computational step is to simply put (I,J) into G. Any 
cycle that forms (a cycle will form if and only if I and J are in the same tree) can always be 
profitably destroyed by removing some other edge (Ki,Li) determined as in Theorem 1. 
(I,J) would be (ii,ji) in that theorem, and z would be positive because the reduced cost for 
(I,J) is negative. Let G^ be G with (I,J) added and, if necessary, (Ki,Li) deleted. Let x^+ 
be X if (Ki,Li) is not deleted, or otherwise the solution x’ of Theorem 1. x^'*' is better than 
X and feasible, but is not in general the G^ basic solution. The latter solution (call it x^) is 
better than x^‘*‘ because it is optimal in P2(G^) while xl'*' is merely feasible, x^ is therefore 
basic and better than x. However, x^ may not be feasible in PI. 

If x^ has one or more negative components, let x2+=x(l-a)-t-xla, where ote[0,l] is 
selected so that x^+ is the closest nonnegative point to xf Let (K 2 ,L 2 ) be an edge in G^ 
whose allocation is driven to 0 at x^+. Since f() is convex and x^ is better than x, x^+ is 
also better than x. x^'*' is also feasible in P2(G^), where G^ is the basis obtained by 
deleting edge (K 2 ,L 2 ) from Gf Let x^ be optimal in G^, so that x^ is better than x^+. If 
x^ also has negative components, continue to define x3+,G^, x^, x^+, etc., until finally x*^ 
is nonnegative, as must eventually happen because an edge is deleted from the basis at each 
step, is then basic, feasible (in PI), and (barring degeneracy) better than x. Since there 
are only finitely many bases, none of which can repeat because the objective function 
decreases at each step, the optimum solution must eventually be found in finitely many 
steps. 

The above algorithm is sufficiently similar to Zangwill's (1970) method of Manifold 
suboptimization that the name has been retained, but there is a difference. The difference is 
that in Zangwill’s method a G-basic, feasible solution that is not optimal would lead to next 
considering a modified basis G'=G-t-(I,J), where (I,J) is some edge at which the reduced 
cost is negative. Unfortunately G' may have a cycle, in which case P2(G') may not have a 
finite optimal solution. The method described above works because every manifold 
considered has a basis that is a forest of trees. 

Manifold suboptimization essentially operates by increasing xu so much that the next 
basic solution may not be feasible, and then fixing the problem. If only feasible solutions 
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are to be considered, then a test for the amount that xij can be increased without causing an 
infeasibility must be found. This is the subject of the next section. 

4. The Extended Simplex Method 

In this section it will be necessary to consider situations where nonbasic variables are 
temporarily positive, as in the Simplex Method. Assume that G is a feasible basis, but that 
the corresponding basic feasible solution is not optimal because KTl is not satisfied for 

some (I,J). Then the reduced cost ?ii-|j.jE(I,J) is negative, so xjj should be increased as in 
Manifold suboptimization. Since the generality will be needed later, suppose that the 
nonbasic variable xjj is fixed at some level b between 0 and bj , and that the balance 
equation (6) is satisfied for each component T of G by z=l when xu=b. If the nonbasic 

variable xij is increased from b to b+A, let the solution of (6) for z be Z(T,A). If f() is a 
good function, Z(T,A) is a continuous, positive function of A with Z(T,0)=1. The 
corresponding multipliers (^(A),|i(A)) are the same as (X*,|i*) except that multipliers for 
nodes in T are multiplied by Z(T,A). Let y(A)=g(|i(A)), let x(A) be the unique solution of 
KT4-KT5 where only basic variables are nonzero except that nonbasic xij=b+A, let A' be 
the largest A such that x(A)>0, and let xkl(^')= 0. In other words, basic variable (K,L) is 

the one first driven to 0 by increasing A. If x(A)>0 for all A>0, let A' and (K,L) be 
undefined. If (K,L) is defined and if its deletion from the basis breaks any cycle that 
introducing (I,J) might form, then the next basis simply replaces (K,L) with (I,J). (K,L) is 
not necessarily defined or part of a cycle even if it is defined, but Algorithm A below 
accounts for these possibilities. The input to the algorithm is a feasible basis G^, 
associated multipliers (X0,|i^), and a nonbasic edge (I,J) with xij=0 and negative reduced 
cost. The output is a different feasible basis G' at least as good as G^, together with 

associated multipliers (X',|i’)- The basic idea of the algorithm is to increase xjj in steps 
until finally the reduced cost is zero. If I and J are in the same tree, xjj is increased in step 
A4 until some edge is deleted. Eventually this operation will put I and J into different trees. 

When I and J are in different trees, either an edge is deleted from one of the trees (A’<A in 
step A5) or else (I,J) is finally added to the basis in the terminal step A9. Here is 
Algorithm A. 

Al)Letk=0, XO=0 

A2) Let X*=XK G=Gk, and b=X^. 

A3) Let Ti be the tree in G^' that includes I and Tj be the tree that includes J. If Ti^Tj 
go to A5. 

A4) Let T=Ti. Determine A' by solving (6) for Z(T,A) and proceeding as in the 
paragraph above. Let (K,L) be the basic edge in T whose variable is first driven to 0, let 
Gk+1=G>^- (K,L), (Xk+l,pk+l)=(X(A'),|i(A')), and Xk+l=Xk+A’. Go to A2. 

A5) Let T=Tp Determine A' by solving (6) for Z(T,A) and proceeding as in the 
paragraph above. Let (K,L) be the basic edge in T whose variable is first driven to 0. Let 

(Ki,Li)=(K,L) and Ai=A'. 

A6) Let T=Tj. Determine A’ by solving (6) for Z(T,A) and proceeding as in the 
paragraph above. If (K,L) is undefined let Aj=oo; otherwise, let (Kj,Lj)=(K,L) and 
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Aj=A'. 

A7) Let A* be the smallest nonnegative solution of the equation 
A.*iZ(Ti,A)=|l*jZ(Tj,A)E(I,J), with A*=o« if there is no such solution. 

A8) Let A'=min{Ai,Aj,A*}. Let ?iik+l=?LikZ(Ti,A’) for ieTi, or ?.ik+l=?iikz(Tj,A') 
for ieTj, or otherwise Define similarly. Let X*^+^=X*^A’. 

A9) If A'=A*, set G’=Gk+(I,J) and (?i', 4 ’)=(>.k+l,^k+l). stop. 

AlO) If A’=Ai, let - (Ki,Li) and k=k+l. Go to A2. 

All) If A’=Aj, let Gk+UG*'- (Kj,Lj) and k=k+l. Go to A2. 

Theorem 3: Algorithm A stops after finitely many steps. The output is a feasible basis 
different from the input, with associated objective function at least as good. 

Proof: A' and (K,L) are well defined in A4 and A5 because row I is included in Tj and 

therefore A' cannot exceed bj. If I and J are included in the same subtree, then the reduced 

cost is negative when A=0 and proportional to Z(T,A), therefore negative for all A. If I 

and J are in different trees, then the reduced cost is negative for A^*, since the functions 

in A7 are continuous in A. Since A’^* in A8, the reduced cost is negative throughout and 
there is no possibility that the objective function might increase. It is also clear that G’ is 
different from G, since G' includes (I,J) and G does not. The only question is whether A9 
is encountered after finitely many steps. 

When an edge is deleted, the effect on the basis is always that there is one more tree 
component and one less edge. Since A2 can only be revisited after deleting an edge, the 
number of edges in G® is an upper bound on the number of times A2 is encountered. Since 
A2 is encountered as frequently as any other step, the proof is complete. QED 

The pivoting operation described above differs from the comparable Linear 
Programming operation in that the number of edges that have to leave the basis in order to 
get (I,J) in is not always 1; it may be any nonnegative integer. It differs essentially in this 
respect from Zan^ill's (1970) Convex Simplex Method (CSM). In the CSM the number 
of basic variables is constant, and every pivot replaces a basic variable with a nonbasic 
variable. When a local minimum is encountered in CSM, the corresponding variable 
remains nonbasic, but positive, and consequently there is no unique solution associated 
with a given basis. In the method discussed above, the number of basic variables is not 
constant, but the idea that nonbasic variables are zero is preserved, along with the one-to- 
one relationship between bases and solutions. This relationship implies that the number of 
pivots is bounded by the number of bases, whereas the number of pivots required by the 
CSM may be unbounded. 

5 . Implementation in the Exponential Case 

Two FORTRAN programs have been written that implement the methods of sections 3 
and 4 in the Exponential case: MANIFO in the case of Manifold suboptimization and 
SIMPLX for the Extended Simplex Method. Each of these programs exploits the forest 
structure of the basis to minimize storage and facilitate computations. The forest is actually 
represented as a single tree by introducing an "earth" node to which all of the components 
of the basis are connected through a root node in each tree; however, in this exposition the 
term "tree" will continue to mean a component of the basis as heretofore. Most of the 
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operations associated with pivoting are similar to those encountered in transshipment 
problems. Since these are well discussed in Bradley, Brown, and Graves (1977), only the 
aspects associated with the nonlinear nature of the objective function will be discussed* 
here. 

Consider first the problem of deleting an edge from the basis. If the edge is in tree T, 
then deleting it amounts to pruning off a branch from T. A reduced version of T will 
remain, along with a new tree T consisting of the pruned off branch. Although equation 
(1) will hold in all trees before the pruning, it generally will not hold in either T or T 

afterwards. The multipliers must be scaled to make (1) hold in T and T', bearing in 

mind that y depends on p. and that nodes in T will have a different z-factor than nodes in 
T. These z-factors can be computed using the formula given in Example 1 of Section 2. 

The computation in T requires the sum pj. This number could be computed by 

jeC(T’) 

summing pj over the appropriate set, but in fact both MANIFO and SIMPLX utilize an 

additional array SDM(k) to store the sum of pj over the set consisting of k and all of its 
column successors in whatever (rooted) tree k belongs to. Thus, if (1,J) is deleted from a 
tree with root node K, and if (say) J is the predecessor of I, then the required sum in T is 
SDM(I) and the required sum in T is SDM(K)-SDM(I). Calculating the z-factor then 

requires a single exponentiation in each tree. After updating (i^,p,y), the allocations in T 
and T' are computed as described in pan b) of Theorem 2.; a preorder traversal array 
permits this to be done in one pass. The allocations themselves are stored in array X, with 
X(k) represenring the allocation from node k to its predecessor if k is a row, or from the 
predecessor to k if k is a column. All arrays are of length n-t-m-t-1 in both programs (node 
length); in fact, the only function requiring mn operations is that of determining the entering 

edge by searching for the largest ratio pjE(i,J)y9q • 

The addition of edge (I,J) to the basis requires the joining of one tree to another, 
possibly after or in conjunction with another pruning operation. If T containing J is to be 
joined to T, T' is first "rehung" so that the root node is J. Next T' is scaled so that 

A.i=pjE(I,J), which requires a logarithm in updating the potentials in T. The unified tree is 
then scaled as in the previous paragraph. It is during the subsequent calculations of X that 
the possibility of negative allocations must be provided for in MANIFO. 

SIMPLX uses all of the arrays in MANIFO plus one additional. The additional array 

W is associated with determining how big A can be before some basic variable becomes 0, 
as will be explained shortly. In the Exponential case the function Z(T,A) satisfies 
ln(Z(T,A))=D(T)A, where 

D(T)=[Xi*5(T,I)-pj*E(l,J)5(T,J)]/ X P j- (7) 

jeC(T) 

Note that D(T) is easy to compute because the denominator is already stored, and also that 
D(T) will be negative in step A5, positive in A6, and negative (since the reduced cost is 

negative) in A4. Multiplying |ij by Z(T,A) is equivalent to subtracting ln(Z(T,A)) from yj, 

so yj(A)=yj(0)-D(T)A for jeC(T). Since this is a linear function of A, allocation xij(A) will 

also be a linear function of A for all edges in T. The array element W(k) contains the rate at 

which the allocation associated with node k decreases with A. The calculation of W is 
equivalent in difficulty to the calculation of X in MANIFO. Given W, it is simple to obtain 
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the limiting value A’ by comparing XAV ratios. The X array can be updated once A' is 
known. 

Both procedures require a basis for a starting point. This is obtained by first including 
all edges of the form (k,k). If there are more columns than rows, this will leave some trees 
consisting of singleton columns. If there are more rows than columns, then row n+k is 
assigned to column k for l<k<n, etc., until finally each row has been included in some 
edge. The resulting basis will be feasible, but unfortunately some edge (i,j) for which 
E(i,j) =0 may be included. The simplest (and the operative) remedy is to simply change 
E(i,j) to some small positive number, since in any case either procedure will quickly delete 
such an edge. 

There is potential for roundoff errors to accumulate in the repeated updating of floating 
point arrays, so NIANIFO and SIMPLX both use double precision arithmetic. Neither 
program incorporates any protection against degeneracy, but so far there have been no 
instances of cycling. 

To give the reader a better idea of the actual operation of the algorithms. Figure 1 
shows the sequence of bases considered by MANIFO in the process of solving a problem 
where m=3 and n=4. The basis is shown as a single tree with row nodes being 
distinguished from column nodes (as they are in MANIFO) by representing row nodes 
with negative numbers. Node 0 is "earth". With the convention of including all edges 
incident to node 0, a basis always has exactly n+m edges. Although X is stored as a node 
length array in MANIFO, the allocations are shown in figure 1 in the more familiar matrix 
form. There are (2^-l)^=3375 bases for this problem, of which MANIFO considers 8 
feasible and 2 nonfeasible before finding one that is optimal. In the course of doing this 
MANIFO requires 18 logarithms or exponentials. These statistics are the same for 
SIMPLX in this small problem; differences emerge only in problems where MANIFO 
introduces more than 1 negative allocation at a time. 

6. Computational Comparisons 

MANIFO and SIMPLX were compared on randomly generated problems to determine 
which program is fastest. Problems generated were such that 
E(i,j) is uniform [0,1] 
bj is uniform [0,10] 

V(j) is uniform [0,100] 

Eleven (m,n) pairs were tested: (10,10), (7,20), (20,7), (20,20), (30,30), 40,40), 
(20,40), (40,20), (50,50), (60,60), and (70,70). Each of eleven (E,b) configurations was 
combined with each of five V configurations to generate a problem. For each of these 55 
problems, the number of nonlinear operations (logarithms and exponentials) NL was 
recorded along with the CPU time T in seconds required to set up the initial guess and 
solve the problem, m, n, NL, and T were then transformed by taking logarithms and a 
linear regression performed. 

There turns out to be very little difference between MANIFO and SIMPLX. The fitted 
formulas are in each case approximately 

T = (.344 sec)(m/46)^-^^(n/122)-^ exp(m/46+n/122) (R2=99%) (8) 

NL=1.4mVn (R2 = 98%) (9) 

All runs were made on the Naval Postgraduate School’s IBM3033AP. Since 
GAMSAIINOS is also available on that system, in one instance a problem with m=20 and 
n=20 was solved using GAMS/MINOS. The time required was 1.88 seconds ("MINOS 
TIME" in the output), which is 47 times as large as the time required by MANIFO (.04 
seconds) for the same problem. Further comparisons were not made because MINOS is a 
general purpose solver that does not exploit the network nature of the constraints. A better 
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comparison would be with the GENOS network optimizer (Mulvey and zfenios (1987)). 
Attempts to accomplish this are in progress, but so far the lack of an exponential function in 
GENOS 1.0 has proved to be a roadblock. 

The reason for keeping track of the number of exponentials and logarithms NL in 
running MANIFO/SIMPLX is that it was initially anticipated that such nonlinear operations 
would require relatively large amounts of time. This turns out not to be the case. When 
(m,n) = (70,70), fomiulas (8) and (9) predict T = .28 seconds and NL = 820. Since an 
IBM3033AP requires only about 4)j,sec for a double precision logarithm or exponential, the 
amount of time spent in nonlinettr operations is orly about (820)(4xl0'^) = .0033 seconds; 
MANIFO and SIMPLX are each mainly occupied with manipulating and scaling the 
various arrays required to represent the current basic solution as a tree, rather than solving 
the balance equation. This is encouraging for applications where the objective function, 
although "gocxJ," might require a numeric rather than analytic solution of the balance 
equation. 
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Solution of a 3x4 Problem Using MANIFO (Continued) 
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Solution of a 3x4 Problem Using MANIFO 
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